Vortex glass transitions in disordered three-dimensional XY models: Simulations for 

several different sets of parameters 
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The anisotropic frustrated 3D XY model with strong disorder in the coupling constants is studied 
as a model of a disordered superconductor in an applied magnetic field. Simulations with the 
exchange Monte Carlo method are performed for frustrations / = 1/5 and / = 1/4, corresponding 
to two different values of the magnetic field along the z direction. The anisotropy is also varied. 
The determination of the helicity modulus from twist histograms is discussed in some detail and 
the helicity modulus is used in finite size scaling analyses of the vortex glass transition. The general 
picture is that the behavior in [Phys. Rev. Lett. 91, 077002 (2003)] is confirmed. For strong (e.g. 
isotropic) coupling in the z direction the helicity modulus fails to scale and it is argued that this is 
due to a too small effective randomness of such systems for the accessible system sizes. 



I. INTRODUCTION 

An applied magnetic field in a type-II superconductor 
will give rise to vortex lines that penetrate the sample. 
A current applied perpendicular to these vortex lines will 
give rise to a force perpendicular to both the current and 
the magnetic field. In a pure system there is nothing 
that hinders the motion of the vortex lines and their mo- 
tion leads to flux-flow resistivity and therefore a loss of 
superconductivity. The presence of point disorder could 
mean a substantial reduction of the mobility of the vor- 
tex lines, but the resistivity would, in the conventional 
picture, nevertheless always be non-zero. 

A vortex glass phase is an alternative possibility 
that was suggested to restore the true superconducting 
statei^ The idea is that the finite disorder strength to- 
gether with the vortex line interaction leads to diverging 
energy barriers against the vortex motion, and thereby 
a vanishing resistivity. This was suggested to take place 
through a continuous transition with universal exponents 
and certain scaling properties. Experimental results in 
support of this picture have been reported^^ but the 
conclusion of a vortex glass phase has also often been 
questioned^ 

There has also been much work on simulations of vor- 
tex glass models. The simplest three-dimensional (3D) 
vortex glass model, that was also the first to be stud- 
ied, is the 3D gauge glass model that includes the dis- 
order through a random vector potential added to the 
phase difference of the superconducting order parame- 
ter. Already the early simulations 7 ^ found strong ev- 
idence for a transition, and with the exchange Monte 
Carlo (MC) technique^ it has been possible both to give 
more convincing evidence for a transition and to de- 
termine the value of the correlation length exponent to 
v= 1.39 ±0.043 

A problem with using the 3D gauge glass as a model of 
a disordered superconductor in an applied magnetic field, 
is the generally recognized fact that the model lacks some 
of the properties and symmetries of the physical system. 
The applied field both breaks the spatial symmetry of 



the system and introduces an additional length scale. In 
a model that properly includes these features one would 
e.g. have the possibility of anisotropic scaling, i.e. differ- 
ent divergences of the correlations parallel and perpen- 
dicular to the applied field. 

Several attempts have recently been done to simu- 
late systems with the correct symmetry. The first pub- 
lished results are from simulations of a frustrated 3D XY 
model with filling / = 1/4 and disorder in the coupling 
constants^ The correlation length exponent was there 
determined to be v = 2.2 even though the quality of the 
data did not allow for any firm conclusions. In a second 
paper by the same author the open boundary conditions 
employed in the first study were changed to standard pe- 
riodic boundary conditions^ The data now rather sug- 
gested v w 1.1, but some quantities still failed to provide 
good scaling. Some aspects of these simulations are dis- 
cussed in Sec. I VII 

Simulations have also been performed with vortex lines 
instead of the phase variables of the XY models A simu- 
lation study of such a vortex line model with strong point 
disorder gives the value v = 0.7, indistinguishable from 
the 3D XY exponent. In that study the pinning energy 
was quite strong which presumably means that most pla- 
quettes are either always occupied or always empty. One 
possible reason for the 3D XY-like exponents could then 
be that the model supports vortex loop excitations (as 
in the 3D XY model) against a background of frozen-in 
field lines from the applied fields 

The present paper is a sequel of Ref. which gave 
the first numerical support for 3D gauge glass exponents 
in vortex glass simulations with the correct symmetry. 
The approach was there to study an anisotropic model 
with much weaker couplings in the field direction than 
in the directions perpendicular to the field. This was 
a natural choice due to experiences from the first order 
transition between the Abrikosov lattice and the vortex 
line liquid. In these simulations 16,17 it has been found 
that the correct behavior required a great flexibility of the 
field induced vortex lines, which could be obtained cither 
with a very large size of the system along the direction 
of the applied field or with weaker couplings between the 
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phase angles in the same direction. As we will see below 
the choice of an anisotropic model turns out to be crucial 
for obtaining convincing scaling collapses. 

Another recent study of the vortex glass transition has 
been done on a model that extends the elastic description 
of a vortex lattice to include dislocations^ The correla- 
tion length exponent of the transition was found to be 
v « 1.3, which within reasonable error bars also is con- 
sistent with 3D gauge glass universality. 

In the present paper we present detailed analyses of 
the frustrated 3D XY model with strong disorder in the 
coupling constants. The paper is an extension of Ref. ITU 
in two respects: (i) The determination of the helicity 
modulus from simulations with twist fluctuations as well 
as analyses of the thermalization and the exchange steps 
in the Monte Carlo simulations are described in consid- 
erably more detail, (ii) Simulations and analyses have 
been done for several different sets of parameters. 

The organization of the paper is as follows: In Sec. 
II we discuss the determination of the helicity modulus 
from twist histograms. Section III deals with the vortex 
glass model and the different sets of parameters used in 
the simulations. In Sec. IV the simulation methods are 
discussed with emphasis on some aspects of the exchange 
Monte Carlo technique, and Sec. V gives the simulation 
results. Section VI, finally, contains a discussion together 
with a short summary. 



II. DETERMINATION OF THE HELICITY 
MODULUS 

The XY model is defined by the Hamiltonian 

(ij) 

where a common choice for the spin interaction is U {(j)) = 
— Jcos(</>). The helicity modulus, which is the stan- 
dard probe of phase coherence in XY models, is defined 
through the response to an applied twist. One way to 
define the twist is to generalize the standard periodic 
boundary conditions #(l,o,o) = #(o.o,o) to 

#(L,0,0) = #(0,0,0) + 

and similarly in the other directions. Here A x is the 
phase mismatch or the total twist in the x direction. One 
may alternatively think about the twist as being spread 
out across the whole system and introduce the twist per 
link, dp = Ap/Lp. The Hamiltonian may then be written 



which gives the correlation function^, 
T 
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With this correlation function the determination of the 
helicity modulus is done in simulations performed with 
zero twist. Note that the derivative in Eq. J2Jl is evalu- 
ated at the minimum of the free energy which typically 
is A = 0. However, in some disordered models there is 
nothing that guarantees that the minimum of the free en- 
ergy is at zero twist. The approach taken here is to study 
such systems with simulations that include the twist fluc- 
tuations as additional dynamical variables. 



A. Twist fluctuations 

There is a well-known duality relation between an XY 
model in the Villain representation and a gas of inter- 
acting charges. In two dimension this is a Coulomb 
gas with logarithmic interactions and in three dimen- 
sions a gas of interacting loops. As observed by several 
authoro 20 ! 21 ! 22 ! 23 the XY model that is dual to a Coulomb 
gas with periodic boundary conditions also includes twist 
fluctuations. Physically, the twist fluctuations are neces- 
sary for the process when a pair of vortices separate, cross 
the boundary and recombine. In the absence of twist fluc- 
tuations such a process gives a configuration where the 
phase rotates by 2tt across the system in the direction 
perpendicular to the vortex separation, as illustrated in 
Fig. ^ The effect is that recombinations of vortices ef- 
fectively is prohibited. Figure |2] illustrates the vortex 
separation in the presence of twist fluctuations in the y 
direction. 

f t t t t t t t t t t t t r t t 
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FIG. 1: The separation of a vortex pair in a system with 
periodic boundary conditions gives a configuration with the 
phase rotating by 2n. 



The helicity modulus is defined through the change in 
the free energy F^A^), or the free energy per site / = 
F/V, as 
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FIG. 2: The separation of a vortex pair in a system with 
fluctuating twist boundary conditions. The twist variable is 
here applied between the top and the bottom rows of spins 
that are connected through the boundary conditions. The 
four panels are for A y — 0, n, 0.85 x 2n, and 2tv, respectively. 
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B. Basic relations 

An alternative means to obtain T M is by first deter- 
mining the free energy, P(A M ). The simulations are then 
performed with fluctuating twists in the /i direction and 
periodic boundary conditions in the other two directions, 

The free energy is obtained from the histogram P(A M ) 
through 

F(A jU )=-TlnP(A jU ), 

and the helicity modulus may be determined from a fit 
of the free energy in a narrow range r of A around zero. 
Dropping the index /i we write, 

F(A) = F Q + ixA 2 , |A|<r. (3) 

This is trivial in principle, but some complications arise 
when this is applied to simulation data with limited ac- 
curacy. The following sections will discuss this question 
in some detail. 



C. Range of A 

Since T is defined as a derivative of the free energy, the 
range of A used for the fit to Eq. @ should be chosen as 
small as possible. To check for the dependence of T on 
the range r we made use of a twist histogram P(A) for an 
ordinary 3D XY model, with L = 8 and T — 2.2 close to 
T c . It is then found that there is a strong dependence on 
r which to a good approximation is T(r) — T(0) ~ — r 2 , 
due to the presence of a A 4 term in F(A). From T(r) 
for small r an extrapolation to r = gives T = 0.1389(3) 
in excellent agreement with the more precise value, T = 
0.13899(8) obtained with Eq. © from a MC simulation 
with the Wolff cluster algorithm. 

D. Disordered systems: unknown A 

We have so far only been concerned with models with 
the known minimizing twist A = 0. The presence of 
disorder may however mean that the minimizing twist 
becomes different for different disorder realizations and 
is not known at the outset, and this turns out to add an 
unexpected complication to the analysis. 

In the case when A is unknown the analysis consists 
of two steps: (i) take some data from a certain range 
around the maximum of P(A) (the minimum of the free 
energy) (ii) fit the free energy from this data to a second 
order polynomial in A to determine T. For this second 
step Eq. has to be changed to 

F(A) =F + ^T(A- A ) 2 , |A-A°|<r. (4) 



When used on simulation data, where statistical fluctu- 
ations are always present, this method happens to give 
values of the helicity modulus that are biased towards 
too large values. 

To illustrate this fact we have again made use of twist 
histograms for the 3D XY model. Even though the min- 
imizing twist is still zero we now take A to be a free 
variable in the analysis. For the complete run which con- 
sists of about 7000 bins there is no discernible effect of 
the randomness, but by constructing twist histograms 
from r aver (say 2-40) consecutive bins the effect becomes 
significant and may be systematically examined. The 
bin size is 2 18 = 262144 sweeps across the system. Fig- 
ure|21shows T(r avcr ) versus l/r av0 r- The values of these 
run-length dependent values T(r ave r) are based on close 
to 7000 different twist histograms constructed from r avor 
consecutive bins: 

P (l) (A;r avcr ) = ^P(A; i + r). (5) 

7~aver * 

T—l 

The message from this figure is that there is a bias in 
determinations of T that are based on too short runs. It 
is also clear that there is a l/r avor -dependence and the 
data may be extrapolated to T(T avC r — ► oo) » 0.1381. 
Since this data is obtained with a finite range, r = 0.0625, 
this value should be compared to, and agrees well with, 
the corresponding value obtained in Sec. Ill CI 
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FIG. 3: When T is determined from Eq. {£) with A as a 
free parameter, the obtained T becomes biased towards too 
large values. This bias decays with l/r aV er where r avor is the 
number of bins used for collecting data, cf. Eq. (|SJ . The data 
is obtained from an isotropic lattice with L = 8 and T = 2.2. 

A clue to the origin of this bias is given by examining 
A which is the location of the minimum of the free en- 
ergy. Since the ground state for the pure 3D XY model 
is a state with zero twist, A = 0, the deviations from 
zero are due to the statistical fluctuations in the twist 
histograms. We find ((A ) 2 ) ~ l/r avcr which is the same 
as the behavior of an average x ol N independent values 
Xt from a distribution with zero mean: (x 2 ) ~ 1/N. 

The key observation is now that T as a measure of the 
curvature of the free energy is inversely related to the 
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width of the distribution. For N samples a;, the width 
of the distribution is characterized by the variance and 
from elementary texts in statistics it is well known that 
CT naivc = lEtf 1 ' "^) 2 gives a biased estimate which 
may be corrected by 



naive ' 

The well known reason for this correction is the use of 
x instead of the true average of the distribution. In the 
analysis of the histogram data the location of A is a 
similar source of error and it is natural to expect the 
same kind of effect in the analysis of the histogram data. 
With T(r aV er) and T(oo) inversely related to c^ aivc and 
a 2 , respectively, and the number of independent samples 
given by N = r avcr /6, we obtain 

T(r avcr ) = — T(oo). (6) 

1 - o/r aver 

Here b is a constant with dimension of time. The expres- 
sion above may also be written 

■YV \ — -YV T _ ^ /''"aver j (7) 

T(T a ver) T(OO) 

and for small values of b/r avel Eq. © becomes 

T(T BV er) = T(00)+67w, (8) 

which explains the rectilinear behavior of T in Fig. 01 T° 
determine the unbiased quantity T(oo) we need to obtain 
T(Taver) for a few values of r avcr and fit that data to one 
of the equations above. 

E. Twist fluctuations in several directions 



"helicity modulus" determined from P x is not the same as 
the proper helicity modulus from the fluctuation formula 
or from P^\ However, from the universality hypothesis 
one would expect scaling of all kinds of quantities based 
on the free energy, and the similar behavior of and 
TW in Fig.Elsuggests that that actually is the case. Here 
we use the standard scaling assumption, 

LT = MtL 1 '"), (9) 

with the reduced temperature t = (T/T c — 1). 




FIG. 4: Comparison of the twist histograms from simulations 
of the pure 3D XY model with twist fluctuations in one and 
three directions. The solid line is the distribution of A x in 
simulations with A y = A z = 0. The dashed line is the same 
quantity obtained with fluctuations in all three directions. 
The data is obtained from a cubic isotropic lattice with L = 8 
and T = 2.2. 



We have now discussed the use of twist fluctuations in 
a single direction and ordinary PBC in the two other. 
The simulations of Ref. were however done in a some- 
what different way with twist fluctuations in all three 
directions. For this discussion we introduce the general- 
ization P (3) (A X , A y , A z ). With the phase angles of the 
XY model discretized to 256 different values, the com- 
puter memory needed to store such histograms rapidly 
becomes enormous. The collected histograms were there- 
fore instead 

^(A.) = ]T ]TP< 3 ) (A,, A„A,), 

Ay A z 

and the analogous P y (A y ) and P Z (A Z ). For the quantity 
defined earlier with twist fluctuations in one dimension 
only we write, 

P^(A x ) = P^(A x ,0,0). 

Figure H shows P^(A X ) together with P X (A X ). These 
two curves are very different and it becomes clear that the 
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FIG. 5: Scaling collapse of the helicity modulii obtained in 
the pure 3D XY model. The upper symbols (solid) are the 
proper helicity modulus, T^' obtained from Eq. © and the 
low symbols are for obtained from P X (A X ). The good 
collapse of the latter quantity confirms the expectation that it 
equally well may be used for examining the critical properties. 
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F. Use in vortex glass simulations 

As discussed above there are some complications in 
the determination of the hclicity modulus from twist his- 
tograms. However, when our interest is only to determine 
the critical properties of a model, two of the above dis- 
cussed complications may be disregarded. If the scaling 
hypothesis is phrased such that the properties of the free 
energy is a function only of the combination tL 1 ^ it is 
clear that the precise method to examine these properties 
is not important as long as it is the same for all system 
sizes. Among other things this means that the choice of 
the range r of P„(A„) used for determining the helicity 
modulus is immaterial. Similarly, the difference between 
the proper helicity modulus and the quantity obtained 
from need not bother us either. The crucial point 
that has to be taken care of is the elimination of the bias 
of Sec. Ill Dl since this bias (as shown in Fig. EI) is different 
for different L. 

When considering disordered systems there is one more 
point that should be taken under consideration. The 
parameter b in Eq. © has the dimension of time and 
may be interpreted as the time between two independent 
measurements. In a disordered system one expects the 
characteristic time to be different for different disorder 
realizations and one would need an average of a num- 
ber of functions with different time constants. However, 
since the correction is linear in b, c.f. Eq. ©, such an 
average has the same functional form, but now with b as 
an average characteristic time. 



III. THE VORTEX GLASS MODEL 

The model we simulate is given by the 
Hamiltoniani^ 

H = - ]T J iM cos(^-^ +A -A iM + ^-«), (10) 

bonds i^i 

where 9i is the phase of the superconducting wave func- 
tion at site i with position of a periodic L x x L y x L z 
lattice, and the sum is over all bonds in directions fj, = x, 
y, z. The size in the x and y directions are the same; 
L x = L y = L. An applied magnetic field in the z direc- 
tion is obtained through the quenched vector potential 
with the choice Ai X — Imfyi, and Ai y = A; 2 = 0. The 
simulations are performed with fluctuating twist bound- 
ary conditions^ which in the duality relation corresponds 
to a vortex line model with periodic boundary conditions. 
We make use of L M twist variables 6^'® in each direc- 
tion and the total twists in the respective directions are 
Afj, = 8$ These variables are updated with the 
usual Metropolis method. We have run simulations for 
four different sets of parameters, summarized in Tab. [I] 
The disorder is put in the coupling constants which are 



chosen as 

Jin = J±(l + £in), fx = x,y, 

Jiz = J\\ (1 + Sip), only sets C and D. 

For set A the £i M are independent variables from a Gaus- 
sian distribution with (e^) = and p = -J (e^) 2 = 0.40. 
For sets B through D were instead from a uniform 
distribution between —1 and 1. Another difference (as 
indicated in the Table) is that the disorder for sets A 
and B is only put on the couplings in the x and y direc- 
tions whereas the sets C and D are also disordered along 
z. The reason for this choice is to facilitate a direct com- 
parison with the simulations in Ref. 



Set 


/ 


J\\/J± 


disorder 


directions 


A 


1/5 


1/40 


Gaussian 


x, y 


B 


1/5 


1/10 


rectangular 


x, y 


C 


1/4 


1/10 


rectangular 


x, y, z 


D 


1/4 


1 


rectangular 


x, y, z 



TABLE I: Four different parameter sets have been simulated. 
Information about the runs are given in 

Tab. El 



IV. SIMULATION METHODS 
A. The exchange steps 

The exchange MC method — also called parallel 
tempering — is an elegant method that makes it possi- 
ble to calculate the correct statistical averages in disor- 
dered systems where the usual MC methods would only 
be stuck in a local minimum. The idea is to simulate 
many different configuration in parallel and, beside the 
ordinary Metropolis MC steps, let the configurations per- 
form a kind of constrained random walk in temperature 
space. These occasional changes in temperature means 
that the configurations sometimes are at higher temper- 
atures where the energy barriers between various local 
minima are low and easily may be overcome. 

Our simulations were done with Nt temperatures, T 
through 7jv t _ i, chosen according to 

/ rp \ m/N T 

T m =T min , m = Q,...,N T -l. (11) 

\ ^ min / 

The values of Nt, T m { ni and T max as well as the number 
of disorder realizations and the length of the runs are 
detailed in Tab. |H| 



B. Check for equilibration 

In spite of its beauty the exchange MC method does 
not alleviate the need for thermalizing the system and 
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Data set 


L 


N d 


N T 


Tmi n 


imin 




Tmax 


A 


10 


600 


12 


0.09 


0.24 


l 


16 


L z /L = 3/5 


15 


600 


24 


0.09 


0.24 


4 


16 




20 


600 


36 


0.09 


0.24 


11 


32 




25 


200 


36 


0.115 


0.24 


17 


48 


A 


10 


900 


12 


0.09 


0.24 


2 


13 


L z /L = 2/5 


15 


900 


24 


0.09 


0.24 


5 


17 




20 


460 


36 


0.09 


0.24 


12 


33 


B 


10 


500 


12 


0.18 


0.40 


3 


13 




15 


500 


24 


0.18 


0.40 


5 


21 




20 


300 


36 


0.18 


0.40 


7 


21 


C 


8 


400 


8 


0.16 


0.38 


1 


12 




12 


700 


16 


0.16 


0.38 


2 


15 




16 


400 


24 


0.16 


0.38 


3 


17 


D 


8 


400 


8 


0.55 


1.10 


1 


13 




12 


600 


16 


0.55 


1.10 


2 


15 




16 


600 


24 


0.55 


1.10 


3 


17 




20 


400 


32 


0.55 


1.10 


4 
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TABLE II: Parameters describing the simulations. For sys- 
tems of size L x L x L z we simulated Nd disorder configura- 
tions with Nt temperatures in the range T m i n < T < T max , cf. 
Eq. |(TTJ. Of the bins corresponding to 2 18 = 262144 sweeps, 
T eq are first discarded and the remaining r K 
for calculating averages. 



Tec are used 



deviate from the true histogram of a hypothetical run of 
infinite length; the difference between a single histogram 
and the true one would give values that are roughly a fac- 
tor of two smaller. The values nevertheless signal large 
fluctuations and conveys a message about a complicated 
phase space with many different local minimas. 
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FIG. 6: The quantity Q M (r) is a disorder averaged mea- 
sure of the difference between the histograms P M (A M ,r) and 
P li (A ll , T max ). The initial decrease down to a constant level 
shows the thermalization of the collection of Nt configura- 
tions. The decrease as r — » r max is there because of correla- 
tions between P M (A M , r) for consecutive r. This data is for set 
A with L z /L = 3/5; the shown quantity is Q — h(Qx + Qy)- 



it is therefore necessary to in some way monitor the ap- 
proach to equilibrium. Since our main quantities from 
the simulations are the histograms P /1 (A i1 ,t) we use 
these quantities in the analysis of the approach to equi- 
librium. The idea is to quantify the similarity of each 
histogram P /J (A At , r) to the last histogram P /J (A /J , r max ) 
which is assumed to be typical of a thermalized system. 
The disorder averaged histogram difference is then de- 
fined as 



£|P„(A„,r)-P„(A,. 



(12) 



The notation [. . .] denotes the disorder averaging. Q = 
\{Qx + Q y ) is shown in Fig. for T = 0.125 (close to T c ) 
and our four system sizes. The decrease of Q at small 
t down to constant levels is an effect of the thermal- 
ization and one can read off the number of bins needed 
for thermalization (r cq in Tab. |nj) from Fig. HJ The de- 
crease of Q as r — > r max is due to similarities between 
P m (A m ,t) and P M (A M ,T max ) that are present because of 
the slow dynamics of the MC simulations. The equilibra- 
tion time r eq is chosen from the time needed to reach the 
constant value of Q(t) whereas r max , the total length of 
the simulations, was chosen to get enough data for the 
extrapolation shown in Fig. |SJ 

The constant level of Q(t) is, especially for the larger 
systems, surprisingly large. To interpret the data cor- 
rectly it should however be kept in mind that the figure 
shows the difference between two histograms that both 



C. Efficiency of the exchange steps 

A common way to monitor the efficiency of the ex- 
change MC steps is to measure the exchange acceptance. 
This is however only a measure of the local mobility of 
the configurations and doesn't answer the more relevant 
question about the efficiency of the algorithm to move 
configurations across a larger temperature range. To 
keep track of all the exchange steps would mean produc- 
ing an enormous amount of data and is therefore usually 
not desirable. A simple method has therefore been de- 
vised that gives the most relevant information with very 
little overhead. The idea is to, for each configuration, 
keep track of the time since the visit at each given tem- 
perature. To that end each configuration is accompanied 
by a vector of integers, v m , with information about how 
long it was since the temperature T m was last visited by 
that very configuration. 

One way to use that information is to examine the 
vectors v m for all configurations that were at the lowest 
temperature at the end of the run. A measure of the dis- 
order averaged time since the last visit at temperature 
T m is shown in Fig. [7| The same figure also shows the 
results from a simple simulation of an unconstrained ran- 
dom walk with the same properties and acceptance prob- 
ability as in the exchange MC. As seen in the figure the 
difference is an order of magnitude, which indicates that 
conclusions about the efficiency of the exchange steps 
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cannot be safely determined from the acceptance ratio 
alone. The reason for the long times needed for a config- 
uration to travel from the highest to the lowest temper- 
ature is presumably that most high-temperature states 
are far away from the phase space regions typical of the 
lowest temperatures, which means that the configuration 
will usually have to undergo many thorough and time- 
consuming reorganizations before it can reach an energy 
compatible with the lowest temperatures. 



x 10 3 



smaller than the errors due to the limited number of dis- 
order realizations. 




FIG. 7: The quantity t v (T m i n ,T m ) is the time, measured in 
number of sweeps, for a configuration to travel from tem- 
perature T m down to T m i n = To. The open symbols show 
the disorder average of the square root of t v based on 600 
disorder configurations. The solid line is an estimate of the 
same quantity from the exchange acceptance. From the large 
difference it is clear that conclusions about the efficiency of 
the exchange steps to make the configurations travel across 
large temperature regions should not be drawn based on the 
acceptance ratio alone. We plot \ftZ since this quantity is 
proportional to the distance for a simple random walk. 



D. Eliminating the bias 
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FIG. 8: The elimination of the bias in Yx(i~ a ver) is done by 
extrapolating to r avcr — + oo with Eq. J7) . The present data is 
from set B, T = 0.2685 closely above T c . 



V. RESULTS 

In this section we report the results from the analysis 
described above, with a number of different sets of pa- 
rameters, cf. Tab. J] The purpose is to check that the 
proposed behavior is a generic feature and is not limited 
to the parameters of Ref. flH but as we have been do- 
ing simulations with several different sets of parameters 
it has not been possible to achieve very high precision 
in the estimates of the critical exponents. The emphasis 
is therefore rather on checking for scaling that is consis- 
tent with 3D gauge glass universality. Generally speaking 
that picture is confirmed, but the new simulations also 
give information about failure of finite size scaling for 
certain sizes and parameters. 



For the following discussion we introduce a notation for 
the disorder averaged helicity modulii in the transverse 
and the parallel directions, respectively, 

T± = i[T x + T,] av , (13a) 

T|| = [TJ av . (13b) 

The procedure used to determine the disorder averaged 
helicity modulus consists of three steps: (i) Determine 
T M (r ave r) for each disorder configuration and several val- 
ues of T avcr by fitting histogram P M (A At ; r aV er) based on 
Tavcr consecutive bins, P /J (A /J ,r) to Eq. ®. (ii) Cal- 
culate the disorder averaged quantities Tj_(r aver ) and 
T || (r aV er), cf. Eqs. JT3J). (hi) Fit this data to Eq. Q 
to obtain the unbiased estimates Tj_ = Tj_(oo) and 
T|| = T||(oo). The last step is illustrated in Fig. |5] 
The error bars are the statistical errors associated with 
the disorder average for each size. It seems that the er- 
rors associated with the extrapolation to zero 1 / r ave r are 



A. High anisotropy, J\\/J± = 1/40 

The results in Ref. were obtained with a rather 
high anisotropy, Ju/Jx — 1/40 and the aspect ratio 
L z /L = 3/5. We have now also performed simulations 
with a smaller aspect ratio, L z /L = 2/5, which is a good 
consistency test since not only the critical exponents but 
also the critical temperature should be independent of 
the aspect ratio. We have also performed additional sim- 
ulations with several other aspect ratios to determine the 
anisotropy exponent. 



1. Varying the aspect ratio 

Figure[5]shows the helicity modulii for the same param- 
eters as in Ref. [FoI but with the aspect ratio L z /L = 2/5. 
We find a nice crossing for LT± at the expected value 
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FIG. 9: Helicity modulii from simulations with aspect ratio 
L z /L = 2/5. The data for LT± in panel (a) all cross at T c w 
0.123 (shown by the arrow) in agreement with Ref . HrL Panel 
(b) for LT|| on the other hand shows the expected crossing 
only for the two largest sizes. The deviation from the scaling 
behavior is another example of the well-known fact that finite 
size scaling often fails for very small system sizes. 



LT 



0.123*i£ The results for the perpendicular quantity 
II also agree with this behavior for the two larger sizes 



but the data for the smallest system, 10 x 10 x 4, is sig- 
nificantly off. This is in line with the general expectation 
that the scaling only should work for rather large system 
sizes. However, somewhat unexpectedly, the scaling in 
LT± prevails even though it fails in the direction paral- 
lel to the applied field. 

Scaling collapses for the two aspect ratios L z /L — 2/5 
and 3/5 are shown in Fig. ^| When discarding T|| for 
10 x 10 x 4 the collapses with T c = 0.123 and v = 1.5 
(from Ref. Il5^ are excellent for both quantities. Note the 
similar shapes of the scaling functions for the two aspect 
ratios. For panel (b) this requires the use of L z instead 
of L on the x axis. Also note that the dependency on the 
aspect ratio is the opposite for L Z T\\ compared to LT±. 



2. The anisotropy exponent 

The above data is consistent with isotropic scaling, the 
anisotropy exponent £ = 1, but to estimate the error 
bars we need a determination of the exponent. The idea 



FIG. 10: Collapse of LT± and L Z T\\ with v = 1.5 and T c = 
0.123 from Ref. Ha for two different aspect ratios, L z /L = 2/5 
and 3/5. The reduced temperature is t = (T/T c — 1). For Ty 
the data for the smallest size, L — 10, has been omitted since 
it appears to be too small to scale, cf. panel (b) in Fig. [5] 



behind finite size scaling is that certain quantities only 
should depend on the fraction £/£, and to generalize this 
concept to anisotropic scaling one has to allow for the 
possibility of two different correlation lengths, £ and £ z , 
that grow in different ways as T c is approached, £ z ~ 
£ c . To do finite size scaling one needs sizes such that 
t;/L oc £z/L z and with the above relation between £ and 
£ z we need to determine the behavior of systems with 
L z ex L/> . For general values of C this gives non-integral 
L z and the common practice is to obtain the appropriate 
data through interpolation of data for neighboring L z - 
values. To make that possible we have thus simulated 
with several different L z : for L = 10 we have used L z = 
5, 6, and 7, and for L = 15 simulations have been done 
with L z = 8, 9, and 10. 

To determine limits on £ the most straightforward test 
would be to repeat the scaling analysis with different val- 
ues of £ and check how the quality of the scaling collapse 
depends on £. Because of the statistical errors in the raw 
data that is however not a very useful technique. A more 
sensitive test is obtained by combining results from anal- 
yses of both and Ty. To do that we focus on how 
the crossing temperatures of L^T± and L 2 ~^T|| depend 
on £. To make the test clean and simple we only make 
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use of two sizes at the time. Fig. ^2 shows the depen- 
dency of the crossing temperatures on ( for sizes L = 15 
and 25. The two different crossing temperatures coincide 
at T s» 0.12 and C ~ 1- Note that the two quantities 
have the opposite dependency on £. This is the key to 
this more precise determination of £, and together with 
a rough error estimate Fig. II II gives £ = 1 ± 0.1. 




1.2 



FIG. 11: The figure shows how the crossing temperature of 
Z/^Tx and L 2 ~^Y| with L — 15 and L = 25 depend on the 
assumed value of f. One set of data is for (L,L Z ) = (25, 15) 
and the other for (L,L Z ) = (15,9 ■ (15/25) c_1 ). The second 
set is obtained by interpolating the results from simulations 
with L z = 8, 9, and 10. Since a crossing of the data in 
both directions should occur at T c the correct value of ( is 
obtained at the crossing of these two sets of data points. This 
gives £ = 1 ± 0.1, strongly suggestive of isotropic scaling. 



way is an artifact of the exchange Monte Carlo method 
since the exchange steps have the effect to give correla- 
tions between results at neighboring temperatures. 
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FIG. 12: Data collapse for / = 1/5 and J\\/Jx = 1/10. The 
collapse is done with the same T c = 0.24 and v = 1.6 for 
both data sets to demonstrate that both quantities may be 
collapsed with the same parameters. 



B. Less anisotropic, J\\/J± ~ 1/10 

The simulations discussed in the previous section are 
for a rather strong anisotropy, J\\/J± = 1/40. It is gen- 
erally expected that the critical behavior should be inde- 
pendent of details as the anisotropy, and we now check 
this expectation with simulations for J\\/J± = 1/10; data 
sets B and C in Tab.[I] Figure IT2"1 shows scaling collapses 
of the helicity modulii for data set B. Beside the weaker 
anisotropy the simulations also differ in that the disorder 
Ei^ is stronger and is now chosen from a uniform rect- 
angular distribution between —1 and 1, corresponding to 

p = J(s%j) = l/\/3 ~ 0.577. In a fit with T c and v as 
adjustable parameters a collapse of LT± gives T c = 0.239 
and v — 1.56 whereas a collapse of LTii gives T c = 0.241 
and v — 1.97. The different values of v is an indica- 
tion of the rather low precision in these determinations. 
In Fig. El we show that it is possible to collapse both 
sets of data with the same parameters, T c = 0.24 and 
v = 1.6. The collapse of LT± is very nice whereas the 
collapse of LTn, especially in a region around T c , is some- 
what worse. However, considering the statistical errors, 
we believe this to be just a statistical fluctuation. The 
fact that several points around T c all deviate in the same 



We have also simulated the same model but with filling 
factor / = 1/4; data set C in Tab.||| The collapse which 
is found in Fig. ^|is excellent and we obtain v = 1.35 
and v = 1.48 from the scaling collapses of LT± and LY\\, 
respectively. 



C. Isotropic system 

The values for v given above, obtained from simula- 
tions with different values of anisotropy and filling fac- 
tor, are within reasonable error bars consistent with 3D 
gauge glass universality, v w 1.39. This seems to rule out 
the possibility that the nice scaling in Ref. was only 
a coincidence. Still, the results presented in the present 
section show that scaling fails when the analysis is ap- 
plied to an isotropic model. This finding is of some im- 
portance since isotropic couplings have been used in sev- 
eral investigationsii*i2ii£ of vortex glass models. These 
papers reach differing conclusions and yet other investi- 
gations fail to find acceptable finite size scaling (private 
discussion). We believe that an understanding of the 
problem to scale our data from isotropic couplings may 
shed light on problems in these other investigations. 
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FIG. 13: Data collapse of LT ± and LT\\ for / = 1/4 and 
J\\/J± = 0.1. The parameters used in the data collapses are 
v = 1.4 and T c = 0.225. 



1. Failure to scale the data 

For simulations of an isotropic system we use the same 
parameters as Kawamura in Ref. [13 but the analysis dif- 
fers from theirs in that we focus on the behavior of the 
hclicity modulii instead of the rms-current. 

Figure I14f a1 shows LT± for the four system sizes 
L = 8, 12, 16, and 20. The data for LT^ weakly sug- 
gests the possibility of scaling and panel (b) shows the 
attempted scaling collapse with T c ss 0.63 and v w 1.50. 
Even though the value of v is in good agreement with our 
earlier findings, the poor quality of the collapse makes it 
impossible to draw any more definite conclusions. Turn- 
ing to LTii/Jii shown in Fig. 1151 we find that it is impos- 
sible to collapse the data since the crossing points for two 
successive system sizes shift systematically to lower tem- 
peratures for increasing L. Beside the failure to scale the 
data it should be noted that LTii/Jii for the isotropic 
case is exceptionally large. For all the other cases we 
had LT||/J|| as 1.0 at T c , but in the isotropic model, this 
quantity is considerably larger in the temperature region 
of interest. 



2. The reason for the failure to scale 

The behavior of T|| in the isotropic system is thus 
clearly different from the anisotropic systems with 
J\\/J± = 1/40 or 1/10. We will now argue that this 
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FIG. 14: Raw data and attempted data collapse of LT± for 
the isotropic system with / = 1/4. The parameters in the 
data collapse in panel (b) are v = 1.50 and T c = 0.63. The 
value of the exponent is consistent with 3D gauge glass uni- 
versality, but the quality of the collapse is not satisfactory. 

is because the disorder in the coupling constants is not 
effective in fully disordering the system for the accessible 
system sizes. 

As a probe of the loss of order we use A° which is 
the position of the minimum of the free energy, ^(A^). 
This quantity has been used before as a measure of the 
effective strength of the disorder^ The disorder fixed 
point was there characterized by (|A°|) = ir/2 which 
corresponds to a uniform distribution between —tt and 
7t. Figure [TBI shows histograms of A° and A° from our 
data and it is clear that the histograms are very different 
from a uniform distribution. Especially the histograms 
of A° are very narrow with |A°/7r| < 0.1 for almost 99% 
of the disorder realizations. For A^ the distributions are 
considerably wider but are still clearly peaked around 
zero. In both cases there is some finite size dependence, 
with a wider distribution for larger system sizes. For 
comparison we also show the corresponding histograms 
for the anisotropic model with J\\/J± = 1/10 in Fig. IT7I 
For the anisotropic case the histograms of A° z are close 
to a uniform distribution; only the data for L = 8 have 
somewhat more weight around zero. This shows that the 
data that exhibits good scaling are from strongly disor- 
dered systems. In contrast, the isotropic model appears 
to be far from the disorder fixed point and we believe 
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FIG. 15: Raw data of LTm for the isotropic system with 
/ = 1/4. Since the crossing temperature for two successive 
sizes shift systematically with increasing L it is impossible to 
collapse the data. The large magnitude of LTy gives addi- 
tional evidence that the isotropic system is different from the 
anisotropic ones. 



that this is at the root of the failure to find a convincing 
data collapse. 

To discuss the physical meaning of we return to 
Fig. [21 which illustrates the relation between the size of a 
vortex pair and the value of the twist variable in the di- 
rection perpendicular to the separation. As the pair sep- 
arates in the x direction the twist Aj, gradually increases. 
At zero temperature the twist is to a good approxima- 
tion proportional to the distance, d, between the vortices, 
A. y = 2ird/L. For the more general situation with sev- 
eral vortices the vortex separation generalizes to the total 
dipole moment of the system of vortices, p x — J2i x i1iy 
where i enumerates the vortices, Xi is the x- coordinate 
of vortex i, and qi is the vorticity (charge). At non-zero 
T the distribution of A y at constant p x will be wider; the 
relevant expressions are given in Ref. Hj, Fot the three- 
dimensional case the dipole moment generalizes to the 
projection of the vortex loops on a certain plane, C xy £& 
The corresponding relation is then A z = 2irC xy / L 2 . 

In a pure system the twist histogram will always be 
symmetric around zero, A° = 0, but the effect of the 
disorder is to favor certain vortex loops between the lay- 
ers and suppress others. The net effect may be a non-zero 
C xy and accordingly a shift of A° z away from zero. Our 
interpretation of the results in Fig. ^j] is therefore that 
the disorder is not strong enough to introduce loops be- 
tween the layers. Note that field-induced vortex lines 
that have a non-vanishing projection on the x-y plane 
also contribute to C xy . The absence of large disorder- 
introduced vortex loops between the layers (or the equiv- 
alent deflection of the field-induced vortex lines) means 
that A° is always close to zero. 

The analyses above suggests that a strong coupling in 
the field direction has the effect to reduce the amount of 
disorder-induced vortex loops between the planes. The 
effect is to get A close to zero which means that the ef- 
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FIG. 16: Histograms of and A° for the isotropic system 
with / = 1/4. In a fully disordered system one expects A° 
to be uniformly distributed between — ir and n, but the figure 
shows that that is not the case for Jm = J± . The peak around 
zero is strongest for A° in panel (a) but is also very clear for 
A® in panel (b). The histograms are calculated on the basis 
of data for all the simulated temperatures. 




FIG. 17: The figures show histograms of A° and A° for 
anisotropy J\\/J± = 1/10 and / = 1/4 in panel (a) and panel 
(b), respectively. The results for the larger sizes, L — 12 and 
16 (circles and squares), are consistent with a uniform dis- 
tribution whereas the distributions for L = 8 (crosses) have 
somewhat more weight around zero. However, it seems that 
such small deviations from perfect disorder (a flat histogram) 
have no discernible effects on the scaling shown in Fig. 1131 
The histograms are calculated on the basis of data for all the 
simulated temperatures. 
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fective disorder is small in the system and we believe that 
this is responsible for the failure of the hclicity modulii 
to scale. Considering the broadening of the histograms 
with increasing L in Fig. 1161 we expect this to be a finite 
size effect, but are presently unable to estimate the size 
where scaling could be expected to set in. 



VI. DISCUSSION 

The use of an anisotropic model in the study of critical 
phenomena with finite size scaling deserves some com- 
ments. To get data with high precision for finite size 
scaling from Monte Carlo simulations, the correlation 
volume should ideally have the same shape as the sim- 
ulation cell. In an isotropic model a cubic simulation 
cell is therefore the best choice and in the general case 
one wants a common value of the fraction £ M /L M in all 
directions. For the model in the present paper with a 
symmetry breaking field there is nothing that guarantees 
that isotropic couplings are best. It is however possi- 
ble to extract some information about the correlations 
from the helicity modulii. With T M as the measure of 
the phase coherence, a larger implies stronger corre- 
lations and thereby a larger correlation length in the \x 
direction. By comparing data for the isotropic model in 
Figs. 1141 and 1151 the fact that T|| is considerably larger 
than Tj_ leads us to conclude that the correlations are 
considerably stronger in the field direction compared to 
the perpendicular direction. One way to reach the goal 
of a simulation cell with the same shape as the correla- 
tion volume would then be to increase the aspect ratio 
L z /L, but a different and more efficient way is to instead 
decrease Jy , the coupling strength in the field direction. 

To get a better understanding of the effect of 
anisotropic couplings on the helicity modulii we have 
made some additional simulations on the ordinary 3D 
XY model (zero field and no disorder) with J z /J = 1/4. 
Since one expects £ M oc y/~J^ the aspect ratio was then 
chosen to be L z /L — 1/2 which gives a simulation cell 
with the same shape as the correlation volume. With 
this value of the aspect ratio the simulations give T/J = 
T z /J z at T c to a good approximation. As shown in 
Fig. E| the same relation holds to a good approxima- 
tion at and close to T c in the simulations of the vor- 
tex glass with J\\/J±_ = 1/10 and L z /L = 1. This 
suggests that the correlations in the different directions 
are about equally strong when the anisotropy is set to 
J\\/J± = 1/10 and that this value therefore is close to 
optimal for the anisotropy in the vortex glass simulations 
with f = 1/5. 

Even though it thus seems that our model is best ex- 
amined with a rather large anisotropy we now turn to 
the results obtained with isotropic couplings. These sim- 
ulations were performed with the parameters of Ref. ^3 
to make it possible to directly compare the results. It 
is however clear that the results are significantly differ- 
ent. Whereas our LT± almost collapse at T = 0.63 with 



v w 1.5 their corresponding quantity, It, collapses for 
the three largest sizes at T g = 0.81 with v = 1.0. Es- 
pecially the different values of the critical temperature 
points to a systematic difference. 

We believe that the reason for this difference is their 
calculation of I rms as the derivative of F(A) evaluated 
at A = 0, rather than at random values of A. From our 
direct determinations of F(A I1 ) we have found that the 
typical structure of this quantity (obtained with param- 
eter set D) is a single minimum of the free energy with a 
shape that in most cases to a very good approximation is 
parabolic, F(AJ = const + T Al (A A1 - A°) 2 /2, where both 
T M and A^ depend on the disorder realization. When the 
derivative is evaluated at A M = one gets 1^ = — T M A°. 
This means that the obtained rms-current is not only a 
measure of the amount of structure in F(A^) but also 
depends on the location of this structure. Against that 
backgound the size-dependence of the distribution of A° 
shown in Fig. 1161 is problematic and we believe this to be 
the reason for the different critical behavior in Ref. ^2 
compared to the results in this paper. This undesired fi- 
nite size will affect the determination of the rms-current 
and we therefore believe that the scaling behavior of It 
in Ref. [l2j is only accidental. 

The failure of the helicity modulus to scale in the 
isotropic model is a related but different question. As 
discussed in Sec. IV C 21 a message from Fig. ^j] is that 
the isotropic model is not sufficiently disordered for the 
simulated system sizes and it seems possible that this re- 
maining order destroys the transition. The broadening 
of the histograms in Fig. 1161 with increasing system size 
would lead to a flat distribution in the limit of large L and 
one would then expect scaling with the 3D gauge glass 
exponents. However, considering the slow widening of 
the histograms as L increases, scaling would presumably 
only be seen for very large systems. 

It should finally be noted that the problem with deter- 
mining the rms-current from the derivative at A = is 
not present in the 3D gauge glass model. The reason is 
that the randomness there is put into the vector poten- 
tial and that a random A then may be absorbed in the 
similarly random Aij . There is then no need to make use 
of a random A and the standard way to evaluate / rms at 
A = is acceptable. 

To summarize: the main conclusion of the present in- 
vestigation is that the vortex glass model is in the same 
universality class as the 3D gauge glass model. This is 
a confirmation of the behavior found in Ref. 0. Still, it 
is found that simulations with isotropic couplings do not 
give any convincing scaling collapse and we argue that 
the reason is that the effective randomness for the acces- 
sible system sizes is too small to give the correct behavior 
of the vortex glass transition. 
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